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Abstract 


In  earlier  papers  the  authors  introduced  and  analyzed  the  calculation 
of  reliable,  a  posteriori  error  estimates  for  finite  element  solutions  and 
discussed  the  design  and  effectivity  of  adaptive  p'ocedures  based  upon 
them.  While  most  of  this  work  concerned  linear  problems,  this  paper  is 
intended  to  show  that  the  same  approaches  remain  also  highly  effective 
in  the  nonlinear  case.  As  a  model  problem  a  planar,  elastic  rod  is  con¬ 
sidered  involving  both  material  as  well  as  geometrical  nonlinearities;  but, 
as  far  as  possible  the  discussion  is  kept  general.  For  nonlinear  problems 
of  this  type  interest  centers  on  an  analysis  of  the  shape  and  features  of 
the  equilibrium  surface.  After  some  general  remarks  about  such  surfaces 
and  the  related  computational  problems,  a  general  continuation  process  is 
summarized  and  some  of  its  extensions  for  determining  limit  points  and 
tracing  stability  paths  are  discussed.  Then  estimators  are  introduced 
for  the  error  along  the  solution  paths  and  of  the  computed  critical  values. 
Experimental  results  for  the  model  problem  show  the  effectivity  of  these 
estimators  and  of  the  adaptive  procedures  based  upon  them.  Then  some  aspects 
relating  to  the  occurrence  and  identification  of  spurious  solutions  are 
discussed.  In  order  to  illustrate  desirable  generalizations  of  these  results, 
the  paper  ends  with  an  outline  of  some  results  for  linear  problems  of  the 
type  that  should  be  achievable  also  in  the  nonlinear  case/- 
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Computational  Error  Estimates  and  Adaptive  Processes 
for  Some  i'lonl  inear  Structural  Pi-oblems^^ 


by 

Ivo  Babuska^^  and  Werner  C.  Rheinboldt^^ 


1 .  Introduction 

In  general,  the  objective  of  a  computational  solution  of  a  physical 
problem  is  to  obtain  an  acceptably  accurate  and  reliable  prediction  of  the 
behavior  of  the  physical  phenomena  under  study.  For  this  the  computation 
should  produce  not  only  an  approximate  numerical  solution  but  also  some 
reliable  information  about  its  accuracy.  More  specifically,  this  accuracy 
should  be  assessed  in  comparison  with  the  solution  of  a  reasonably  general 
mathematical  model  of  the  problem  and  in  terms  of  a  specified  error  norm 
that  matches  the  requirements  of  the  study. 

The  effectivity  of  such  an  estimation  may  be  judged  by  the  effectivity 

index 

(11)  0  =  estimated  error 

true  error 

In  practice  it  is  usually  more  important  for  e  to  be  close  to  one  than  that 
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0  2i  1 •  Moreover,  it  is  essential  that  e  converges  to  one  when  the  errors 
converge  to  zero,  so  that  for  an  accuracy  of  say  5S  or  10%  the  value  of 
[e-ll  is  expected  to  be  less  than,  say,  0,1  or  0.2,  The  requirement  of 
designing  error  estimators  with  these  properties  for  realistic  classes  of 
problems  and  various  different  norms  certainly  represents  a  demanding  task. 

In  order  to  achieve  the  stated  objective  of  the  computation  adaptive 
processes  appear  to  be  optimal  tools.  In  fact,  it  is  beginning  to  be  a 
widely  accepted  observation  that  for  realistic  problems  it  is  rarely  feasible 
to  design  numerical  processes  which  reliably  and  effectively  achieve  the 
desired  accuracy  in  the  prescribed  norm  and  yet  which  do  not  utilize  some 
form  of  adaptivity.  Unfortunately,  the  term  "adaptation"  remains  to  be  rather 
ill-defined  and  is  often  used  very  loosely.  A  clarification  of  this  concept 
may  require  close  attention  to  results  in  such  fields  as  automatic  control 
theory,  artificial  intelligence,  learning  processes,  etc.  (see  eg.  [22], 

[33],  [42],  [44]).  In  our  context,  the  availability  of  precisely  defined, 
reliable  error  estimators  appears  to  be  central  to  the  design  of  effective 
adaptive  procedures.  The  effectivity  has  to  be  measured  here  in  terms  of 
the  size  of  the  class  of  problems  for  which  solutions  within  a  specified 
range  of  accuracy  are  obtained  with  minimal  cost. 

These  general  concepts  for  the  numerical  solution  of  certain  practical 
problems  have  been  discussed  by  us  in  several  recent  papers.  More  specifi¬ 
cally,  for  the  finite  element  solution  of  certain  problems  in  structural 
mechanics  we  introduced  and  analyzed  the  calculation  of  reliable,  a  posteriori 
error  estimates,  [4  ],  [ 5  ],  [ 6  ],  [ 9  ],  [13],  as  well  as  the  design  and 
effectivity  of  adaptive  procedures  based  on  these  estimates,  [7  ],  [8  ],  [10], 
[11],  [12],  [46]. 

Up  to  now  essentially  all  this  work  concerned  linear  problems  in  one 
or  two  space  dimension.  Not  unexpectedly,  for  nonlinear  problems  the 
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situation  is  considerably  more  difficult  and  much  remains  to  be  done.  In 
the  second  part  of  [12]  some  results  were  given  for  the  finite  element 
solution  of  a  particular  nonlinear  model  problem  in  which  showed 
that  a  posteriori  error  estimates  and  adaptive  approaches  are  also  highly 
effective  in  such  nonlinear  cases  even  though  the  corresponding  theory  is 
as  yet  far  from  complete. 

This  paper  is  intended  to  be  a  continuation  of  this  work  on  nonlinear 
problems.  Such  nonlinear  problems  show  special  features  not  present  in  the 
linear  case.  In  particular,  they  involve  usually  a  number  of  intrinsic 
parameters  and  interest  centers  not  so  much  on  determining  a  few  specific 
solutions  for  fixed  parameter- values  but  to  assess  the  behavior  of  these 
solutions  under  general  variations  of  the  parameters.  In  structural  analysis 
the  parameters  may  characterize  load  points  and  load  directions,  material 
properties,  geometrical  data,  etc.,  and  the  set  of  all  solutions  depending 
on  these  parameters  has  been  called  the  equilibrium  surface  of  the  structure 
(see  [37]).  This  equilibrium  surface  provides  considerable  insight  into  the 
behavior  of  the  structure  and  its  stability  properties;  we  refer,  for  instance 
to  [30],  [38]  for  further  discussions  and  various  examples.  From  a  numerical 
viewpoint  the  question  then  is  to  analyze  computationally  the  shape  and 
characterize  features  of  this  equilibrium  surface.  For  this  a  principal 
tool  is  a  general  form  of  continuation  process  which  allows  a  trace  of  a  priori 
specified  as  well  as  intrinsically  defined  paths  on  the  surface. 

In  the  analysis  of  a  structural  problem  by  the  finite  element  method 
we  can  compute  only  approximate  points  along  a  path  on  the  solution  surface 
of  some  discretized  form  of  the  original  problem.  Then  —  in  line  with  the 
earlier  stated  objectives  of  a  practical  computation  --  we  are  faced  with  the 
need  for  assessing  and  controlling  the  errors  along  an  entire  segment  of 
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such  a  path.  In  addition,  since  we  are  interested  in  specific  features 
of  the  equilibrium  surface,  such  as,  the  location  of  turning  points, 
bifurcation  points,  stability  boundaries,  etc.,  we  require  also  estimators 
of  the  error  between  such  specific  points  and  their  computed  approximations. 

Some  a  priori  estimates  of  this  type  have  been  given  recently  in  [16],  [17], 
[18].  But,  as  in  the  linear  case,  we  require,  of  course,  computable  and 
reliable  a  posteriori  estimators  which  can  then  form  the  basis  of  adaptive 
procedures  for  computing  with  a  prescribed  accuracy  the  desired  segments 
of  solution  paths  on  the  equilibrium  surface. 

In  the  first  part  of  this  paper  we  consider  a  strongly  nonlinear,  one¬ 
dimensional  rod-problem  and  show  for  it  that  effective  a  posteriori  error 
estimates  and  adaptive  procedures  can  be  constructed  and  successfully 
integrated  into  a  general  continuation  process  for  tracing  paths  on  the 
equilibrium  surface.  In  addition,  we  introduce  error  estimators  for  the 
computed  location  of  turning  points  along  such  paths.  Finally,  various 
numerical  results  show  the  effect! vity  of  all  these  approaches  and  ideas. 

The  results  given  here  are  as  yet  restricted  in  their  applicability 
but  should  be  capable  of  generalization  to  broader  classes  of  problems.  In 
order  to  illustrate  where  such  generalizations  should  lead  us  we  present 
in  Part  II  of  this  paper  an  outlook  on  some  further  results  for  linear 
problems  of  the  type  we  would  hope  to  achieve  ultimately  also  for  the  non¬ 
linear  case.  In  particular,  this  includes  the  design  and  application  of  an 
adaptive  finite  element  solver  for  certain  two-dimensional  elliptic  problems 
and  the  extension  of  the  use  of  error  estimators  and  adaptive  procedures  to  some 
parabolic  problems.  For  further  summaries  of  the  specific  results  presented 
in  this  paper  we  refer  to  the  outlines  included  at  the  beginning  of  each  one 
of  the  two  parts. 
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Part  I:  Nonlinear  Problems 

In  this  part,  consisting  of  sections  2  to  9  of  this  paper,  we  discuss 
various  aspects  of  the  computational  solution  of  a  class  of  nonlinear 
problems  in  structural  mechanics.  More  specifically,  we  consider  as  a 
model  problem  a  planar,  elastic  rod  involving  both  material  as  well  as 
geometrical  nonlinearities.  But,  as  far  as  possible,  the  discussion  is 
kept  general  and  applicable  to  broader  classes  of  problems.  As  indicated 
in  the  introduction,  for  nonlinear  problems  of  this  type  we  are  interested 
in  analyzing  the  shape  and  features  of  the  equilibrium  surface.  Accordingly, 
we  include  here  some  general  remarks  about  such  surfaces  and  about  the  related 
computational  problems  involving  solution  paths,  bifurcation  and  turning 
points,  stability  boundaries,  etc..  Then  a  general  continuation  process  is 
summarized  and  some  of  its  extensions  for  determining  limit  points  and 
tracing  stability  paths  are  discussed.  In  line  with  the  comments  in  the 
introduction  we  introduce  estimators  for  the  error  along  the  solution  paths 
and  of  the  computed  critical  values.  Experimental  results  for  the  model 
problem  show  the  effect! vity  of  these  estimators  and  of  the  adaptive  pro¬ 
cedures  based  upon  them.  Finally  we  discuss  some  aspects  relating  to  the 
occurrence  of  so-called  spurious  solutions  which  occur  frequently  in  dis¬ 
cretized  nonlinear  problems. 
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2.  A  Model  Problem 

As  a  model  problem  we  consider  the  planar  deformations  of  a  non-li nearly 
elastic  rod.  More  specifically,  large  deformations  as  well  as  material  non- 
linearities  are  allowed.  A  comprehensive  treatment  of  the  global  behavior 
of  such  rods  was  given  in  [1].  These  rods  can  suffer  not  only  flexure  but 
also  compression  (or  extension)  and  shear. 

The  basic  configuration  of  the  rod  is  sketched  in  Figure  1.  We  restrict 
ourselves  to  the  case  of  a  symmetric  shape  where  it  suffices  to  consider 
only  one  half  of  the  rod,  say,  the  left  half  in  our  picture.  Then  the  con¬ 
figuration  of  the  rod  is  described  by  two  functions 


(2.1) 


r:  [0,n  q:  [0,1]  -  Ma(s)lU  =  1 


where  r  defines  the  axis  of  the  deformed  rod  and  £  the  direction  of 

the  cross-section  (expressing  the  effect  of  the  shear  stresses).  The  parameter 

s  represents  the  arc-length  of  the  rod-axis  in  the  straight  reference  con- 
1  2 

figuration.  Let  e  ,  e  be  the  global  basis  vectors  as  shown  in  the  figure  and 

2 

u  the  angle  between  £  and  e  .  Then,  with 
(2.2)  U(u)  = 


cos  u 

sin  u 

-sin  u 

cos  u 

we  define  the  local  coordinate  system  by 


(2.3) 


£(s)  =  U(u(s))'’^e\  £(s)  =  U(u(s))^e^ 


and  the  strains  n  by 


(2.4) 


1  +  C(s) 
n(s) 
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r'{s)  =  U(u(s)) 


Let  f(s)  and  m(s)  denote  the  resultant  force  and  moment,  respectively, 
in  the  cross-section  at  the  point  s.  Since  the  load  is  constant,  the 
equilibrium  equations  have  the  form 

(2.5)  f'  =  0,  m’  +  I’  1  =  0 

where  primes  denote  differentiation  with  respect  to  s.  Hence,  we  have 


(2.6) 


and 


(2.7) 


f(s)  =  m(s)  =  M(s)(e^  X  e^) 


0  =  M'(s)  +  (p^ 


r'(s)  =  w(s)  +  [u(u(s))  (;;)f 


1  +  Us)'\ 
n(s)  /  ■ 


We  introduce  the  constitutive  equations 

M(s)  =  a(u'(s)), 

(2.8)  C(s)  =  b((e^)^  U(u(s))  C^)}  =  b(-x  cos  u(s)  +  v  sin  u(s)), 

n(s)  =  c((eb^  U(u(s))  C))  =  c(v  cos  u(s)  +  x  sin  u(s)), 

n 


where  a,  b,  and  c  are  given  real  functions.  With  the  abbreviation 


(2.9) 


b(u,x,v)  =  (x  sin  u  +  \)  cos  u)(l  +  b(-x  cos  u  +  v  sin  u)) 


+  (x  cos  u  -  V  sin  u)  c(x  sin  u  +  v  cos  u) 
our  equation  (2.7)  now  assumes  the  form 

(2.10)  ^  a(u')  +  b(u,X.v)  =  0. 

The  fixed  left  endpoint  and  the  symmetry  requirement  result  in  the  homo¬ 
geneous  boundary  conditions 

(2.11)  u(0)  *  u(l)  =  0. 

Note  that  for  a(t)  =  a^t  and  b(t)  h  c(t)  5  0  the  problem  reduces  to 
the  classical  Bernoulli -Euler  problem 

(2.12)  a^u"  +  X  sin  u  +  v  cos  u  =  0,  u(0)  =  u(l)  =  0. 

Equation  (2.10)  is  strongly  nonlinear  when  a(t)  ^  a^t  and  only  mildly 
nonlinear  when  a(t)  =  a^t  as  in  (2.12).  Mildly  nonlinear  problems  of 
this  form  have  been  studied  by  many  authors.  For  some  recent  results  about 
finite  element  approximations  for  such  equations  see,  for  instance,  [16J,  [17],  [18] 
and  the  references  given  there. 

In  order  to  define  the  specific  constitutive  equations  used  in  our 
numerical  experiments  let 

(itt  for  t  <  0 
gi(t)  =  arctan  -^t,  g2(t)  =  \ 

t  +  2t^  for  t  >  0 


(2.13) 
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Then  the  three  cases  of  (2.8)  considered  in  this  paper  are  as  follows: 

(i)  (Strongly  nonlinear)  a{t)  =  g^(t),  b(t)  =  g2(t),  c(t)  =  g^(t) 

(2.14)  (ii)  (Mildly  nonlinear)  a(t)  =  |-  t,  b(t)  =  g2{t),  c(t)  =  g^(t) 

(iii)  (Classical  Euler  Rod)  a(t)  =  j  b(t)  =  0,  c(t)  =  0 

Note  that  from  (2.4)  and  (2.8)  it  follows  that  the  shape  of  the  axis 
of  the  rod  is  specified  by 
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3.  Remarks  on  Equilibrium  Surfaces 

Our  model  problem  (2.10)/(2.11)  involves  the  two  parameters  x,  v 
representing  the  components  of  the  resultant  force.  For  any  fixed  x,  v 
a  solution  u  =  u(s)  of  the  boundary  value  problem  defines  an  equilibrium 
configuration  of  the  rod.  For  practical  applications  it  is  rarely 
sufficient  to  determine  only  the  equilibria  for  a  few  given  loads.  Instead 
we  are  interested  usually  in  the  changes  of  the  equilibrium  position  under 
specific  changes  of  the  load. 

Frequently  in  structural  mechanics  for  a  fixed  direction  the  load- 
intensity  is  used  as  a  parameter.  For  example,  in  our  case  we  may  consider 
the  family  of  loads  -xe^  acting  in  the  direction  of  the  axis  of  the  rod 
in  its  undeformed  configuration.  If  X  varies  say  from  x  =  0  to  some 
maximal  value  we  obtain  in  this  way  a  curve  of  equilibria  of  the  rod.  Often 
this  curve  is  parametrized  in  terms  of  x.  But  at  turning  points  this 
parametrization  breaks  down  and  hence  it  is  usually  better  to  consider  the 
curve  in  the  form  x  =  X(t),  u  =  u(t)(s)  with  another,  more  suitable 
parameter  t.  Clearly,  for  different  load  directions  or  starting  points  we 
obtain  a  different  curve. 

In  our  case  the  problem  depends  on  the  two  parameters  x  and  v. 
Accordingly,  it  turns  out  that  all  these  curves  form  a  two-dimensional 
surface  in  (x,v,u)-space,  the  so-called  equilibrium  surface  of  the  problem. 
In  general,  this  surface  has  a  complicated  shape  but  its  consideration 
provides  considerable  insight  into  the  equilibrium  behavior  of  our  structure. 

In  order  to  illustrate  these  concepts  it  may  be  helpful  to  consider 
for  a  moment  a  particular  frame-work  model  of  our  rod  (see  also  [30],  p.  291). 
Two  rigid  rods  of  length  1  each  are  flexibly  jointed  as  shown  in  Figure  2 
and  a  linear  spring  with  spring-constant  k  is  trying  to  hold  the  structure 
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in  its  straight  reference  configuration.  For  simplicity,  let  k  =  1/2. 

Any  equilibrium  position  is  characterized  by  the  value  of  the  angle  y, 
and  it  is  easily  seen  that  the  equilibrium  equation  has  the  form 

(3.1)  y  -  X  sin  y  -  1  V  cos  y  =  0. 

Hence,  the  equilibrium  surface  here  is  a  two-dimensional  surface  in  the 
three-dimensional  space  with  coordinates  x,  v,  y.  It  is  a  ruled  surface 
for  which  the  generators  have  constant  values  of  y.  Figure  3  shows  level 
lines  of  a  projection  of  part  of  this  surface  in  the  direction  of  the  y-axis. 

In  addition,  Figures  4  and  5  illustrate  the  cross  sections  of  the  surface  with 
planes  x  =  const  and  v  =  const,  respectively.  The  latter  two  figures  show 
clearly  that  the  point  x=l,  v=0,  y=0  is  special.  More  precisely, 
we  have  here  a  bifurcation  point  where  two  equilibrium  curves  for  the  load 
direction  (1,0)^  are  intersecting.  The  particular  shape  of  the  surface 
near  this  point  represents  a  so-called  cusp-catastrophe.  The  reason  for 
this  name  is  obvious  from  Figure  3. 

Evidently,  the  equilibrium  surface  characterizes  much  better  the  complete 
behavior  of  the  framework  than  the  usual  traces  of  the  equilibrium-curves 
corresponding  to  fixed  load  directions.  Furthermore,  the  surface  may  be 
partitioned  into  subsets  of  stable  and  instable  equilibria.  As  usual,  an 
equilibrium-solution  is  called  (locally)  stable  if  it  represents  a  local 
minimum  of  the  total  energy.  In  our  case,  it  is  readily  seen  that  instability 
occurs  when 

(3.2)  1  -  X  cos  y  +  V  sin  y  <  0. 

In  Figure  3  the  upper  part  of  the  visible  surface  consists  of  stable  equilibria. 


Figure  ^ 
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The  boundary  of  the  stability  region  is  the  envelope  of  the  generators. 

This  stability  boundary  is  shown  as  a  dashed  line  in  Figure  4.  It  con¬ 
sists  of  the  turning  points  with  respect  to  the  load  parameter  \  on  the 
curves  for  constant  values  of  v.  These  points  represent  the  buckling 
points  of  the  structure  for  the  particular  values  of  v.  Thus,  at  one 
glance  the  surface  provides  us  with  information  about  the  stable  equilibria, 
buckling  loads,  etc.  for  our  framework. 

This  indicates  that  for  a  deeper  understanding  of  the  behavior  of  a 
structure  we  should  analyze  the  form  of  the  equilibrium  surface.  Of  course, 
the  choice  of  the  parameters  x  and  v  entering  into  the  definition  of  our 
surface  is  basically  arbitrary  albeit  natural.  Practical  considerations 
show  that  the  load-component  x  is  essential  for  the  study  of  the  buckling 
behavior  of  the  framework.  At  the  same  time  it  is  necessary  to  introduce 
one  or  several  parameters  describing  possible  perturbations.  There  are 
theoretical  reasons  to  consider  a  particular  minimal  number  of  such  per¬ 
turbation  parameters,  but  we  shall  not  enter  into  that  question  here. 

In  view  of  these  observations  the  aim  of  our  analysis  is  to  determine 
the  form  of  the  equilibrium  surface.  In  our  simple  example  this  is  easy 
since  no  discretization  is  involved.  But  in  general,  we  can  compute  only 
points  which  are  approximately  on  the  solution  surface  of  a  discretization 
of  the  original  problem.  In  line  with  the  earlier  stated  requirements  for 
any  practical  computation  this  raises  the  question  of  assessing  the  errors 
conmitted  here. 

Analogous  to  the  framework  the  equilibrium  surface  of  the  rod-problem 
(2.10)  is  a  two-dimensional  surface  in  the  three-dimensional  space  with 
coordinates  x,  v,  u.  Suppose  that  u{x,v)  and  u(x,v)  denote 
the  exact  solution  and  the  computed,  approximate  solution,  respectively,  for 
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given  parameter  values  x,  v.  An  obvious  possibility  for  the  error  assess¬ 
ment  then  is  the  requirement  that  |ju(x,v)  -  u(x,v)||£  e.  However,  this 
is  reasonable  only  if  both  u  and  u  exist  for  the  particular  parameter 
values.  This  is  by  no  means  guaranteed  as  Figure  6  shows  which  provides 
a  sketch  of  possible  curves  u(*,v^)  and  u(*,v^).  For  x  <  Xj  both 
u  and  u  exist  but  for  x  >  x.^.  there  is  no  value  of  u{x,v^)  which 
can  be  compared  with  the  computed  value  u(x,Vq).  Hence,  it  will  be 
necessary  to  consider  different  measures  of  the  error  between  the  exact 
and  computed  solutions.  A  possibility  is  to  ask  for  an  estimate  which 
guarantees  that  for  any  computed  u(x,v)  there  exists  an  exact  solution 
u(x,v)  with 

(3.3)  |X  -  xl'^  +  |v  -  v(^  +  llu(x,v)  -  u(x,v)|]  _<  c 

and  suitable  exponents  a,  S.  Clearly,  this  avoids  the  difficulty  with 
the  earlier  type  of  estimate.  In  addition  to  estimates  such  as  (3.3)  other 
more  specific  error-assessments  are  often  desired.  For  example,  it  is 
frequently  of  considerable  practical  interest  to  measure  the  error  |Xy  -  x.j.1 
in  the  parameter  values  corresponding  to  the  turning  points  of  the  exact 
and  computed  solutions. 

In  connection  with  this  discussion  it  should  be  noted  also  that  the 
equilibrium  surface  of  the  discretized  problem  may  consist  of  more 
connected  components  than  the  surface  of  the  exact  problem.  In  other  words, 
there  may  be  numerical  solutions  which  do  not  approximate  exact  solutions. 

Such  spurious  solutions  have  been  observed  in  many  contexts  (see  eg.  [20],  [28] 
and  the  references  given  there).  From  our  viewpoing  the  question  then  is 
to  identify  the  segments  of  solution  paths  on  the  equilibrium  surface  of  the 
discrete  problem  for  which  the  error  estimates  are  below  an  acceptable  threshold. 
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Figure  6 
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4.  A  Finite  Element  Approximation  of  the  Rod  Problem 

We  return  again  to  the  rod  problem  (2.10)/(2.n ) .  The  finite  element 
solution  of  this  boundary  value  problem  is  based  on  the  standard  variational 
form 


(4.1) 


fl 

[-a{u')v'  +  b{u,x,v)v]ds  =  0 

•'0 


where  v  varies  over  a  given  space  of  test  functions. 

A  mixed  finite  element  method  was  used  to  approximate  the  solutions  of 
(4.1).  More  specifically,  let 


(4.2)  a:  0  =  Sq  <  s.|  <  $2  <  • . .  <  =  1 ,  n  =  n(A) 

be  a  given  partition  of  the  interval  [0,1].  Then  on  the  i-th  interval 
=  [s^_-j,s^.],  1  i  £  n,  of  A  the  solution  is  approximated  by 


(4.3)  u.(s)  =  ^  y.  .  i^.(t.(s)),  V  s  £  I.,  1  1  i  1  n, 

where  , . . .  ,ii,p  are  the  standard  Legendre  polynomials  on  [-1,+1] 

and 


(4.4)  T.(s)  =  (s  -  (s..^  +  s.)). 


h.  =  Si  - 


M-l 


At  the  interior  nodes  of  A  the  functions  and  their  derivatives  up  to  order 
d-1  are  to  match  continuously  (2d  _<  p*l).  Thus,  together  with  the  two 
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boundary  conditions  at  Sq  and  we  have  to  add  d(n-l)  +  2  constraints 
to  (3.1).  In  order  to  take  account  of  the  corresponding  constraints  for 
the  test  functions  an  equal  number  of  Lagrange  multipliers  have  to  be  intro¬ 
duced.  Hence,  altogether,  we  have  (p+l)n  +  d(n-l)  +  2  variables.  In 
the  standard  finite  element  method  the  use  of  the  corresponding  elements 
results  in 

(4.5)  (p+l)n  -  d(n-l)  -  2, 

degrees  of  freedom.  This  number  is  used  later  in  all  estimations  of  the 
total  work. 

The  mixed  method  of  this  form  was  selected  for  our  example  computations 
because  it  allows  for  a  flexible  use  of  elements  of  various  degree  and 
smoothness.  With  this  it  was  possible  to  assess  among  other  points  the 
amount  of  work  for  a  fixed  accuracy  corresponding  to  different  degrees  of 
approximation  and  varying  smoothness.  High  order  elements  are  used  fre¬ 
quently,  for  instance,  in  the  P-version  of  the  finite  element  method,  (see 
eq.  [14],  [15],  [41]). 
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5.  A  Continuation  Method 

In  structural  mechanics  continuation  methods  are  usually  viewed  in  a 
narrow  sense  as  methods  for  following  numerically  a  specific  load  curve 
parametrized  by  a  load  intensity.  As  we  discussed  in  section  3  the  solution 
set  of  structural  problems  generally  constitutes  some  surface  or  manifold. 

For  example,  for  our  rod-problem  (2,10)/(2.11 )  the  equilibrium  surface  is  a 
two-dimensional  surface  in  (.\,v,u)-space.  For  a  thorough  understanding  of 
the  equilibrium  behavior  of  a  non-linear  structure  we  are  interested  in 
determining  the  principal  features  of  the  equilibrium  surface,  as,  for 
example,  the  location  of  critical  boundaries  and  bifurcation  points  where 
stability  is  lost,  etc,. 

For  a  numerical  analysis  of  a  given  equilibrium  surface  we  need  to 
consider  continuation  methods  in  a  broader  sense  as  a  collection  of  numerical 
procedures  for  accomplishing  tasks  of  the  following  type: 

(i)  Follow  numerically  any  curve  on  the  surface  specified 
by  a  particular  combination  of  parameter  values. 

(ii)  Determine  on  such  a  curve  the  critical  points  where 
stability  may  be  lost. 

(5.1) 

(iii)  From  any  such  critical  point  follow  a  path  in  the 
critical  boundary. 

(iv)  On  any  curve  specified  by  (i)  determine  the  location  of 
bifurcation  points  and  the  paths  intersecting  at  that 
point. 


Methods  for  these  and  related  tasks  have  been  developed  by  various 
authors  in  recent  years  (see  eg.  [23],  [24],  [27],  [29],  [31],  [32],  [34],  [39], 
[40]  and  the  references  given  there).  We  shall  not  go  into  details  here  but 
sketch  only  some  of  the  principal  features  of  these  methods  as  they  apply 
in  the  present  setting. 
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Generally,  a  (finite-element)  discretization  of  a  structural  problem 
leads  to  a  finite  dimensional  non-linear  equation  of  the  form 

(5.2)  F(y.p)  =  0 

where  y  e  is  a  vector  of  state  variables  (eg.,  deformations),  p  e  a 
vector  of  parameters,  and  F;  r”  x  r™  a  given  function.  For  ease  of 

discussion  we  shall  assume  here  that  m  =  2.  The  set 

(5.3)  E(F)  =  {(y,p)  e  r"  X  r2  1  F(y,p)  =  0} 

of  all  solutions  of  (5.2)  is  the  equilibrium  surface  of  our  discretized 
problem. 

Let  DF(y,p)  denote  the  derivative  of  F;  that  is,  the  n  x  (n+2)- 
( Jacobi  an)  matrix  of  all  partial  derivatives  of  F,  We  assume  that  the 
so-called  regularity  set 

(5.4)  R(F)  =  {(y,p)  e  r"’  X  R^  I  rank  DF(y,p)  =  n) 

is  not  empty.  For  example,  in  the  case  of  the  framework  of  Figure  2  we  have 
y  e  r\  p  =  (A,v)^  £  R^  and 

(5.5)  DF(y,p)  =  (4k  -  2x  cos  y  +  v  sin  y,  -2  sin  y,  -cos  y). 

The  three  partial  derivatives  in  this  1  x  3  matrix  DF(y,p)  can  never 
be  zero  at  the  same  time.  Hence,  the  rank  of  that  matrix  is  always  one  and 
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3 

the  regularity  set  R(F)  is  the  entire  space  R  . 

On  the  equilibrium  surface  (5.3)  we  may  define  a  curve  by  specifying 

some  relation  between  the  two  components  a,  v  of  the  parameter  vector 
T  2 

P  =  (x,v)  e  R  .  For  example,  suppose  that  we  fix  a  linear  relation  of 
the  form 

(5.6)  q'^(p-p°)  =  qi(A-X°)  +  q2(v-v°)  =  0 


where  q  =  (qi.q2)^  is  a  given  non-zero  vector  and  p°  =  (x°,v°)  a 
specified  vector  of  parameter  values.  Hence,  we  are  now  interested  in  the 
solutions  of  (5.2)  for  which  (5.6)  holds;  in  other  words,  we  are  considering 
the  expanded  equation 


(5.7) 


F(y,p)  = 


=  0. 


In  our  examples  the  parameter  vector  p  represents  a  load.  Hence, 
p°  is  a  reference  load,  (5.6)  defines  a  load  direction  and  the  solutions 
of  (5.7)  are  the  equilibria  of  the  structure  in  this  particular  load 
direction. 

The  derivative  of  F  is  the  (n+1)  x  (n+2)  matrix 


(5.8) 


DF(y,p) 


F{y,p) 


and  it  is  readily  seen  that  the  regularity  set 
of  R(F);  that  is 


R(F)  of 


F 


is  a  subset 
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(5.9) 


R(F)2I  R(F), 


It  is  important  to  note  that  in  general  there  are  points  in  R(F)  which  do 
not  belong  to  R(F).  For  example,  for  the  framework  of  Figure  2  suppose 
that  we  are  interested  in  tracing  a  curve  defined  by  v  =  v°.  Hence,  in 
(5.6)  we  use  q  =  (0,1)^,  p°  =  (0,v°)^,  and  (5.8)  has  the  form 


(5.10) 


DF(y,p)  = 


3 

The  matrix  has  rank  two  everywhere  in  R  except  on  the  straight  lines  de¬ 
fined  by  X  =  (-1)\  y  =  iiT,  i  =  0,+^l,+2,...  and  variable  v.  Thus, 

R(F)  is  no  longer  equal  to  R  as  was  the  case  with  R(F). 

The  intersection  of  the  equilibrium  surface  E(F)  with  the  straight 
lines  where  DF(y,p)  has  only  rank  one  are  the  points  with  the  coordinates 
X  =  (-1)\  V  =  2(-l  )^''’^i7r,  y  =  itr,  i  =  0,+l,+2,...  .  In  the  region 
where  |y|  _<  17/2  the  only  such  point  is  the  bifurcation  point  X  =  1,  v  =  0, 
y  =  0  of  Figure  4.  Note  that  this  point  became  a  singular  point  only  by 

considering  the  family  of  curves  v  =  v°  on  the  equilibrium  surface. 

Different  families  of  curves  generate  different  singular  points.  For  instance, 
if  we  trace  the  family  of  curves  x  =  x°  for  given  x°,  then  in  our  strip 
|yi  ^  77/2  the  only  points  not  in  R(F)  have  the  coordinates  x  =  77/2, 
y  =  +  77/2,  V  =  -2.  This  effect  is  clearly  visible  in  Figures  4  and  5. 

For  a  fixed  vector  q  the  points  of  the  regularity  set  R(F)  are  of 
principal  interest  for  our  continuation  processes,  since,  under  simple 
conditions  about  F,  there  exists  a  unique  solution  curve  of  (5.7)  through 


any  such  point.  More  precisely  we  have  the  following  existence  result  (see 
[32]); 
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n  9  n 

Theorem:  Let  F:  r"  x  r”  be  continuously  differentiable  and  DF  locally 

o  2  ^ 

Lipschitzian.  With  given  q,p  e  R  form  the  expanded  mapping  F  of  (5.7). 

Then  for  any  (y°,p°)  £  R(F)  there  exists  a  unique,  continuously  differentiable 
solution  y:  J  C  R*^,  p:  J  C  R^  -*■  R^  such  that 


R(y(s),p(s))  =0,  Vs  £  J,  y(0)  =  y°,  p(0)  =  p° 

(5.11) 

(y(s),p(s))  e  R(h,  l|(y'(s),p'(s))l  I2  =  1,  VseJ 


where  the  open  interval  J  C  R^  is  maximal  under  set-inclusion.  Moreover, 
if  s  is  a  finite  endpoint  of  J  then  for  s  s,  s  e  J  either  the 
points  (y(s),u(s))  tend  to  the  boundary  of  R(F)  or  | ly(s)  ,u(s))] | ^ 
goes  to  infinity. 

The  effect  of  this  theorem  is  clearly  visible  in  Figure  4.  Through  any 
point  in  the  shown  region,  except  the  bifurcation  point  X  =  1,  u  =  0,  y  =  0, 
there  is  exactly  one  solution  curve  which  terminates  only  at  a  bifurcation 
point.  Note  that  at  x  =  tt/Z,  y  =  it/2  the  v  =  v°  curves  in  R^  do  not 
intersect  although  in  the  figure  this  appears  to  happen. 

For  the  discussion  of  a  continuation  process  for  computing  a  curve  (5.11), 
it  is  useful  to  use  the  abbreviation  x  =  {y,u)  for  the  vectors  (y,u)  e 
It  can  be  shown  that  on  the  regularity  set  R(F)  there  exists  a  unique  tangent 
mapping  T  such  that 


T:  R(F)  DF(x)Tx  =  0,  MTxI]^  =  1, 


(5.12) 


>  0. 
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In  other  words,  for  any  point  x  e  R(F)  the  vector  Tx  of  length  one  is 
tangent  to  the  unique  solution  curve  through  that  point.  The  determinant 
condition  in  (5.12)  was  introduced  to  fix  the  orientation  of  that  tangent 
direction. 

Most  modern  numerical  continuation  processes  for  tracing  a  curve  (5.11) 
are  of  the  predictor-corrector  type.  But  the  methods  differ  in  the  choice 
of  the  parametrization  of  the  curve.  We  shall  discuss  here  a  locally  para¬ 
metrized  process  of  the  form  presented  in  [32],  [34], 

0  1  k 

Let  X  ,x  ,...,x  ,  0.  be  the  approximate  points  along  the  curve 

k+1 

which  have  been  computed  so  far.  The  computation  of  the  next  point  x 
has  schematically  the  following  form: 

k  k 

1 .  At  X  compute  the  tangent  vector  Tx  . 

2.  Choose  the  index  i  ^  ''  '^be  variable  x. 

which  is  to  serve  as  the  current  continuation  parameter. 

k  k 

3.  Choose  a  steplength  tj^^i  along  the  Euler  line  x  +  tTx  . 

ored  k  k 

4.  Compute  the  predicted  point  x^  =  x  +  t.^,  Tx  and  the 
corresponding  "corrector-constant" 

5.  With  as  starting  point  apply  Newton's  method  to  the 

corrector  equation 


Fx 

If  "no  convergence"  then  replace  tj^,^  by  ^  and  go  to  4, 

If  "convergence"  then  use  the  last  iterate  as  the  next  continuation 
point 


For  the  computation  of  the  tangent  vector 


Tx" 


note  that 


Ji 
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X  /  i  T  k  / 

(5.14)  det  I  1  =  [(e  )  Tx  ]  det  [  1  ,  i  =  l,...,n+2 


where  e\...,e'’^^  are  the  natural  basis  vectors  of  Hence,  with 

1  T  k 

any  index  i  such  that  (e  )  Tx  ^0  the  tangent  vector  may  be  computed 
as  follows 


(5.15) 


DF(x'^) 

(e'Y 


.n+2 


Tx  =  a 


a  =  [sign  det 


1  i  T  k 

.•  T  I  ]/[sign  (e  )  Tx^]. 

[{e')  ’ 


For  k^  1  we  use  the  index  i  =  i|^_i  of  the  continuation  variable  during 
the  previous  step  while  for  k  =  0  the  index  of  the  original,  user-provided 
continuation  variable  is  taken. 

In  line  with  (5.14)  the  index  i  =  \  of  the  variable  x^.  which  is 
to  serve  as  the  current  continuation  variable  is  chosen  such  that 


(5.16)  l(e  Tx*^|  =  max  Ke’)"*^  Tx*^]. 

i=l , . . . ,n+2 

This  represents  the  variable  for  which  the  basis  vector  e^  forms  the 

smallest  angle  with  the  tangent  Tx  .  Note  that  —  unless  the  step  from 
k-1  k 

X  to  X  was  unreasonably  large  —  the  choice  of  i  =  in  (5.15) 

should  guarantee  that  the  matrix  indeed  is  nonsingular. 

In  order  to  determine  a  steplength  we  proceed  as  in  the  case  of  ODE- 
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solvers  and  introduce  the  quadratic  polynomial 

(5.16)  q(t)  =  +  tTx*^  +  t^ 

with 

(5.17)  w*"  =  (Tx"  -  -J-  (x'^-x'""'')),  =  ||x^-x''‘''|l2 

"k  k 

for  which 

(5.18)  q(0)  =  x*^,  q'(0)  =  Tx*^,  q{-h|^)  =  x'^"\ 

The  distance  between  and  the  curve  is  to  be  smaller  than  some  tolerance 

>  0«  Asymptotically  for  sufficiently  small  hj^  and  t|^^^  this  re¬ 
quirement  may  be  approximated  by 

(5.19)  -  ^1(^1^112-  Vl 
whence 

/^k+1 

The  polynomial  (5.16)  is  numerically  ill -defined  when  the  angle 


(5.21) 


=  arc  cos  [(Tx*')^  (j|j-  (x*^  -  ^ ))]  e  [O.tt] 
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is  small,  that  is,  when  the  path  is  fairly  straight.  This  makes  it 
advisable  not  to  attempt  to  use  q(t|^^^)  in  place  of  as  the 

predicted  point.  The  norm  of  w  may  be  evaluated  easily  in  the  form 


(5.22)  |lw''ll2  =-^  Isin^l  . 

The  choice  of  the  "corrector-constant"  Y|^  in  line  4  of  (5.13) 

defines  the  point  on  the  path  (5.11)  which  is  to  be  approximated  by  the 

1/ 

corrector  iteration.  More  specifically,  suppose  that  x  is  an  approxi¬ 

mation  of  the  point  x{S|^)  on  the  curve,  then  for  given  and  sufficiently 
small  t  and  As  the  equation 

(5.23)  (e’"''^)^  (x(s^+As)  -  (x^t Tx*") )  =  y^ 


defines  a  one-to-one  correspondence  between  t  and  As.  Evidently,  the 
choice  t  =  AS  underlying  the  inequality  (5.19)  is  a  natural  but  by  no 
means  required  goal.  If  we  replace  x(s)  in  (5.23)  by  q(t)  t^®n 
AS  =  t|^^i  implies  the  use  of  the  corrector  constant 


(5.24) 


.  Vi  [(.''kjT  I 


Obviously,  Y|^  should  be  evaluated  in  double-precision  to  reduce  the  possible 
effect  of  subtractive  cancellation  when  the  curve  is  nearly  straight. 

For  the  definition  of  the  steplength  this  leaves  the  choice  of  the 
tolerance  .20).  A  study  of  some  aspects  of  this  question  has 

been  given  in  [19J.  For  our  discussion  here  it  suffices  to  mention  only  a 


simple  but  rather  effective  technique. 


For  any  real  numbers  a  <  b  we  define  the  piecewise  linear  function 


(5.25) 


u(t;a,b) 


a  if  t  _<  a 

t  if  a_<t_<b 

b  if  t  >  b 


The  steplength  calculation  uses  as  input  the  angle  of  (5.21),  the  pre¬ 
diction  error 


(5.26) 


6k  =  ll(x'''^+tk  Tx'"-'')  -  x"" 


at  the  last  step,  and  the  last  step  hj^  defined  in  (5.17).  With  these 
quantities  the  calculation  is  essentially  based  on  (5.20)  and  proceeds  as 
follows: 


1.  If  then  6:  =  k,  go  to  4  [straight  line  case] 

else  if  then  6:  =  go  to  4 

[hairpin  case] 


(5.27) 


2.  6:=  u(|p;  6^^^)  [adjustment  of  relative  tolerance] 

k 


3.  6^:  =  vj( - - ,-4.<^)  [adjustment  of  0] 

2|sin  la  '  - 


4.  :  •  u(9h,, 


L 


[adjustment  of  step] 
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Here  0  <  <  tt  are  bounds  for  the  angle  ctj^  and  k  is  the 
maximally  allowable  step-increase,  that  is,  we  want  to  ensure  that  the 
step- factor  6  is  bounded  by 


(5.28) 


1 

—  <  0  <  K. 

K  —  — 


In  line  with  this  the  tolerance  bounds  in  line  2.  of  (5.27)  are  defined  by 


(5.29) 


sin  i 
mm  2  2 

K 


max’ 


max 


=  2k 


sin  2 


mm 


The  condition  6„.  <  places  some  restriction  on  the  choice  of  a  .  , 

mm  max  mm 

3nd  K.  Experience  has  shown  that  a  =  0.05,  a  ,  =  tt/2,  and 
max  mm  max 

<  =  3  are  effective  choices  for  a  large  range  of  problems.  The  actual 

step  is  enforced  to  fall  between  the  bounds  0  <  t  .  <  t  .  Here  t  . 

mm  max  mm 

depends  on  the  accuracy  of  the  machine  and  t  reflects  the  desired 

max 

density  of  points  along  the  curve. 

Line  5.  of  (5.13)  requires  some  provisions  for  monitoring  the  con¬ 
vergence  of  the  corrector  iteration.  For  Newton-type  methods  it  has  been 
found  satisfactory  to  declare  non-convergence  if  either  one  of  the  three 
conditions  is  true 


(i)  I  |Gx*'’'^[  I  _>  9|  |Gx*^’'^”^  1 1 ,  for  some  j  1 
(5.30)  (ii)  I  |x'^’'^'-x*^’'^'^  1 1  2  ‘’I  »  for  some  j  _>  2 


k  i 

where  lx  *'’}  denotes  the  sequence  of  corrector  iterates,  (ji  is  a  suit¬ 
able  constant,  say  t  =  1.05,  and  an  integer,  say,  j'  =5.  On 

iTIaX  iTIaX 

k  i  k+1 

the  other  hand,  the  iterate  x  is  accepted  as  the  next  point  x  if 

(5.31)  ( I  iGx*^”^!  I  £  e.|)  and  ( |  lx*^’'^'-x*^’'^'”^  1 1  _<  ^2  ^3!  ) »  '^or  some  j 

with  given  tolerances  £i,£2.e2  which  depend  on  the  machine  and  the  problem. 

As  stated  in  (5.13),  in  the  case  of  non-convergence  the  predictor  step 

is  halved  and  the  corrector  is  restarted  with  the  corresponding  value  (5.24) 

of  the  corrector-constant  Y|^.  Of  course,  here  it  is  required  that 

J  1  ^min  action  is  required.  On  the  other 

k*^!  k  i 

hand,  if  convergence  is  declared  with  x  =  x  as  the  last  iterate 
then  a  control  program  has  to  decide  whether  a  new  predictor/corrector  step 
is  needed  or  not. 

The  algorithm  described  here  has  been  used  extensively  and  with  excellent 
success  on  a  wide  range  of  problems.  In  particular,  it  is  well  suited  for 
the  task  (i)  of  (5.1).  For  practical  use  it  is  desirable  to  incorporate 
special  features  which  allow,  for  instance,  for  the  detection  and  computation 
of  any  target  point  on  the  curve  where  a  specific  variable  has  a  given  value, 
or  for  the  detection  and  computation  of  a  limit  point.  The  latter  feature 
will  be  discussed  briefly  in  the  next  section. 
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6.  Determination  of  Limit  Points  and  Critical  Boundaries 

In  this  section  we  discuss  briefly  some  of  the  algorithmic  aspects 
relating  to  the  tasks  (ii)  and  (iii)  of  (5.1),  that  is,  to  the  determination 
of  critical  points  on  a  continuation  curve  and  the  trace  of  a  path  in  the 
critical  boundary.  As  before  let  (5.2)  be  the  given  system  of  non-linear 
equations  which  defines  the  equilibrium  surface  E(F)  of  (5.3)  for  our  dis¬ 
cretized  problem. 

n  2  1 

Suppose  that  there  exists  a  total  energy-functional  tt:  R”  x  R  R 
T  n  2 

with  Dy7r(y,p)  =  -F(y,p),  V  y  e  R”,  p  e  R'^.  Then  an  equilibrium 
(y°,p°)  e  E(F)  is  stable  if  Tr(y,p°)  has  a  proper  local  minimum  at  y  =  y° 
and  hence  if  DyF(y°,p°)  is  negative  definite.  Let 

(6.1)  C(F)  =  {(y,p)  £  R^  X  rP  I  D^F(y,p)  is  singular); 

then  the  so-called  critical  boundary  C(F)nE(F)  contains  the  (relative) 
boundary  of  the  set  of  stable  equilibria  on  the  equilibrium  surface  E(F), 
although  the  two  sets  are  not  necessarily  equal.  The  tasks  (5.1)  (ii)  and 
(iii)  concern  the  computation  of  points  and  curves  in  the  critical  boundary 
C(F)  n  E(F),  Note  that  this  set  is  also  well-defined  if  there  is  no  total 
energy  functional  tt.  But  then  D^F(y,p)  is  non-symmetric  and  stability 
may  be  lost  at  points  not  in  C(F)n  E(F). 

As  in  the  previous  section  we  consider  curves  in  E(F)n  R(F)  defined 
by  an  augmented  equation  (5.7).  As  we  saw  before  the  regularity  set  R(F) 
of  the  expanded  function  F  is  contained  in  R(F).  More  specifically,  it 
can  be  shown  that 


(6.2) 


R(F)  =  LR(F)\C(F)]  UL(F,q) 
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where 

(6.3)  L(F,q)  =  {(y.p)  e  r"  X  1  dim  ker  D^Fly.p)  =  1,  DpF(y,p)  ^  rge  DyF(y,p) 

In  other  words,  either  a  point  (y,p)  e  R(F)  on  our  curve  does  not  belong 
to  C(F)  or  it  is  in  the  set  L(F,q).  As  the  notation  already  indicates 
this  set  L(F,q)  depends  only  on  F  and  the  vector  q  used  in  (5.6)  to 
specify  the  curve.  Moreover,  one  can  show  that 

(6.4)  L(F,q)  =  {X  e  R(F)  |  (e^)''^Tx  =  0,  j  =  n+l,n+2}. 


In  other  words,  the  points  (y,p)  e  L(F,q)  on  our  curve  are  exactly  those 

points  of  R(F)  where  the  tangent  of  the  curve  is  orthogonal  to  the  2- 

dimensional  subspace  of  the  parameter  variables.  Thus,  in  standard  termi- 

nology>the  points  of  L(F,q)  are  the  limit  points  of  F  in  R(F)  or  more 

precisely  the  limit  points  of  F  with  respect  to  the  direction  q. 

Suppose  now  that  our  continuation  process  of  Section  5  is  applied  to 

trace  the  curve  defined  by  (5.7)  and  a  specific  starting  point 

x°  =  (y°,p°)  e  R(F)  for  which  F  x°  =  0.  As  shown  in  [32]  the  computation 

will  continue  as  long  as  the  curve  remains  in  R(F).  In  order  to  detect 

on  the  curve  a  limit  point  of  F  we  have  to  look  for  points  of  L(F,q). 

t  T 

For  this  it  is  only  necessary  to  monitor  one  component  (e  )  Tx  of  the 
tangent  vector  Tx  along  the  curve,  provided  the  index  t,  n+1  £  _<  n+2, 

was  chosen  such  that  the  vector  q  has  a  non-zero  £-th  component;  that  is. 


(6.5) 
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This  follows  readily  from  the  characterization  of  L(F,q)  and  the  form 
of  K 

If  our  starting  point  x°  e  R(F)  represented  a  stable  equilibrium 
then  the  first  limit  point  encountered  on  our  path  indicates  a  point  where 
stability  may  be  expected  to  be  lost.  For  structural  problems  this  is 
typically  a  point  where  buckling  may  occur.  Let  x  =  x{s)  e  R{F),  s  0, 
x(0)  =  x°,  be  the  segment  of  the  exact  curve  which  we  are  tracing.  Then  the 
first  limit  point  is  the  first  zero  s  =  s*  >_  0  of  the  real  function 

(6.6)  <^(s)  =  Tx(s).  s  >  0. 

where  the  index  I,  n+1  ^  ^ n+2,  was  chosen  such  that  (6.5)  holds. 

If  there  exists  a  total  potential  and  (}>  changes  signs  at  s  =  s* 
then  at  that  point  of  the  curve  one  eigenvalue  of  D^F(x(s*))  becomes 
positive  and  stability  is  certainly  lost.  Limit  points  where  $  changes 
sign  are  called  simple  limit  points.  They  may  be  characterized  in  terms 
of  the  behavior  of  the  second  derivative  D^^F{x{s*))  but  there  is  no  need 
to  discuss  this  here  (see  eg.  [35]).  We  note  also  that  zeroes  of  (j)  where 
the  sign  does  not  change  are  numerically  ill-conditioned  and  cannot  be 
detected  easily.  Accordingly,  in  practice,  interest  centers  almost  exclusively 
on  simple  limit  points. 

The  continuation  process  produces  a  sequence  of  points  x°,x^,... 

approximating  our  curve  and  we  expect  to  encounter  a  first  pair  of  points 
k-1  k 

x  ,  X  where 

sign  (e"^)^  Tx*^’^  f  sign  (e^)^  Tx*^. 


(6.7) 
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This  signals  the  nearness  of  a  simple  limit  point  x*  e  L(F,q)  and 
the  task  (5.1)  (ii)  requires  a  closer  computation  of  x*.  For  this 
several  algorithms  are  available  (see  the  references  given  in  [35]).  A 
very  simple  but  reliable  approach  is  to  apply  a  root-finder  to  the  function 
<t>.  More  precisely,  we  consider  the  secant  line 

(6.8)  z(t)  =  (l-t)x*^"^  +  tx*^,  0  £  t  _<  1 

k-1  k 

between  x  and  x  and  denote  by  i  the  index  of  the  maximal  component 
k  k-1 

in  modulus  of  x  -  x  .By  construction  of  our  continuation  process 
Newton's  method  applied  to  the  system 

/  \ 

(6.9)  \  =  0 

I  (e'y  (x-z(t))j 

should  converge  to  a  point  on  the  curve.  Let  x{t)  denote  the  accepted, 
final  Newton-iterate  corresponding  to  the  starting  point  z(t),  0^  t^  1. 
Then  the  real  function  ({)(t)  =  (e"^)^  Tx(t),  0  _<  t  j<  1,  is  well-defined  and, 
by  (6.7),  it  has  different  signs  for  t  =  0  and  t  =  1.  Hence,  a  standard 
root-finder  can  be  applied  to  4).  In  particular,  experience  has  shown  that 
the  code  given  in  [21]  works  here  very  well. 

In  structural  problems  it  is  rarely  sufficient  to  determine  just  one 
limit  point.  If,  say,  the  structure  is  expected  to  buckle  at  that  point 
then  we  are  interested  in  determining  the  change  of  that  buckling  point 
with  a  variation  of  the  parameters.  The  standard  computational  approach 
for  this  is  to  repeat  the  entire  process,  that  is,  to  define  a  neighboring 
curve  by  changing  the  starting  point  and  the  augmenting  equation  (5.7)  and 
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to  detennine  on  this  new  curve  the  first  simple  limit  point.  Obviously, 
this  may  be  a  very  time-consuming  method.  Moreover,  it  is  usually 
extremely  difficult  to  construct  neighboring  curves  which  lead  to  a 
sequence  of  reasonably  closely  spaced  limit  points. 

This  raises  the  question  about  the  feasibility  of  methods  for  a 
direct  trace  of  a  curve  in  the  critical  boundary  C(F)n  R(F)  through  our 
first  simple  limit  point  x*.  Several  such  methods  have  been  discussed 
in  [35]  and  it  would  lead  too  far  to  go  into  details  here.  Briefly,  one 
of  the  processes  works  with  two  augmented  systems  of  the  form 


(6.10) 


G^Cy.P.P) 


F(y.p) 


(q^)^  (p-p) 


=  0,  ^  =  1,2 


2  2  12 

where  q  =  q  e  R  and  q  e  R  is  orthogonal  to  q.  In  general,  these 
systems  define  two  families  of  curves  on  the  equilibrium  surface.  Then  the 
concept  of  the  process  is  to  start  at  the  first  limit  point  x*  with  a  trace 
of  the  curve  through  x*  defined  by  Gp  After  one  or  at  most  a  few  con¬ 
tinuation  steps  we  switch  to  the  curve  defined  by  G2  and  determine  on  it 
the  first  simple  limit  point  which  represents  the  next  point  on  our  desired 
curve  in  the  critical  boundary  (see  Figure  7).  A  repetitive  application  of 
this  concept  requires  attention  to  some  technical  details  for  which  we  refer 
to  the  cited  article  [35]. 

As  an  example  of  the  application  of  this  process  we  consider  a  clamped, 
thin,  shallow,  circular  arch  which  has  been  used  as  a  test  case  by  various 
authors  (see  eg.  [25],  [26],  [43]).  Let  U  and  W  be  the  radial  and  axial 
displacements,  R  the  arch-radius,  A  the  cross-sectional  area,  h  the 
thickness,  and  E  Young's  modulus.  With  the  dimensionless  displacements 
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u  =  U/n,  w  =  W/h,  the  total  potential  energy  —  non-dimensionl ized  by 

2 

dividing  by  EAR(h/R)  --  is  given  by 


(6.11) 


-6, 


- 

de 


pude , 


-e. 


Here  p  =  p(e)  is  the  dimensionless  radial  load  and  ct^  ,a2  are  dimensionless 
constants.  Each  end  is  assumed  to  be  pinned,  that  is,  we  have  the  boundary 
conditions 


(6.12) 


u(+  e^)  =  0, 


w(+  9q)  =  0, 


d^u 


de 


The  finite  element  approximation  introduced  in  [43]  was  applied.  More 
specifically,  we  used  a  uniform  mesh  with  eight  rlements,  6^  =  15°  and 
the  constants  a-j  =  3.8716  x  10'  ,  02  =  8.2752  x  10"  corresponding  to 
the  data  in  [26].  The  following  four  load  functions  p  =  p(iJ,v)  were  chosen: 


(1)  P(u,v)  = 


u(l-v)  +  4iiv,  for  elements  4  and  5 
u(l-v),  otherwise 


(2)  p(y,v) 

(6.13) 

(3)  p(u,v) 


fu(l-v)  +  4ijv,  for  elements  2  and  3 
u(l-v),  otherwise 

ju(l-v)  +  8iiv,  for  element  4 
|y(l-v),  otherwise 


(4)  p(u,v) 


u(l-v)  +  8uv,  for  element  2 
u(l-v),  otherwise 


Hence,  in  all  case  the  average  load  is  y. 

If  the  boundary  conditions  are  disregarded  there  are  45  displacement 
variables  and  two  control  parameters.  In  order  to  simplify  the  comparison 
for  different  boundary  conditions,  all  variables  were  always  numbered  con¬ 
secutively  from  1  to  47.  Thus,  X21  denotes  the  radial  displacement  of 
the  center  point.  Of  course,  once  the  boundary  conditions  are  taken  into 
account  there  remain  only  39  displacement  variables. 

For  each  load  function  a  basic  curve  on  the  equilibrium  surface  was 
defined  by  fixing  a  value  of  v.  More  specifically,  for  the  symmetric  load 
(6.13)  (1)  the  value  v  =  0  was  used  while  in  all  other  cases  v  =  0.005 
was  taken.  These  curves  were  followed  numerically  until  a  first  limit  point 
was  detected.  From  this  limit  point  the  corresponding  curve  in  the  critical 
boundary  was  traced  by  the  method  sketched  above.  The  results  are  shown 
in  Figure  8.  The  curves  in  the  critical  boundary  for  the  load  functions 
(6.13)  are  labeled  (1)  through  (4).  They  represent  the  location  of  the 
buckling  points  of  the  arch  under  the  chosen  load  regimes.  Note  the  con¬ 
siderable  influence  of  the  location  of  the  load  maxima  for  equal  average 
loads.  The  curves  (2),  (3),  (4)  for  the  asymmetric  loads  emanate  from  the 
bifurcation  point  on  the  basic  curve  defined  by  v  =  0.  On  these  curves 
the  structure  buckles  asymmetrically.  On  the  other  hand,  for  the  symmetric 
load  (1)  the  "critical  curve"  intersects  at  the  limit  point  of  the  curve 
for  V  =  0,  and  represents  points  where  the  structure  undergoes  symmetric 
"snap  through"  buckling.  It  is  easily  seen  that  this  is  an  instable 
situation  which,  in  practice,  is  not  realizable  except  under  very  controlled 
conditions . 

The  example  shows  that  this  approach  can  provide  a  wealth  of  information 
about  the  buckling  behavior  of  a  structure.  Moreover,  it  turns  out  that  the 
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computational  work  for  tracing  one  of  the  curves  in  the  critical  boundary 
is  only  mildly  larger  than  that  for  tracing  one  of  the  basic  curves.  Thus, 
this  type  of  analysis  of  a  structure  is  certa'nly  practically  feasible. 
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7.  A  Posteriori  Error  Estimates  and  Adaptive  Approaches 

In  Section  4  we  introduced  a  particular  finite  element  formulation 
for  our  model  problem  (2.10)/(2.11 ) .  This  formulation  is  rather  flexible 
and  leaves  the  user  considerable  freedom  in  the  choice  of  various  data-items, 
including,  in  particular,  the  mesh  (4.2),  degree  p  of  the  elements  and 
smoothness- index  d  of  the  solution.  Now  we  face  the  following  two  tasks: 

(i)  Construct  a  posteriori  error  estimates  which  --  for 
given  input-data-i terns  —  characterize  accurately  and 
reliably  the  error  caused  by  the  finite-element 
solution  as  measured  in  some  specified  norm. 

(ii)  Design  an  adaptive  procedure  for  selecting  the  controlling 
data-items  defining  the  solution  so  as  to  achieve  an 
acceptable  solution  with  the  prescribed  accuracy. 

These  two  tasks  are  closely  tied  together.  Without  reliable  error 
estimators  it  is  impossible  to  design  effective  adaptive  procedures  for 
solving  the  problem  within  the  prescribed  error-range.  At  the  same  time, 
it  can  be  shown  theoretically  and  practically  that  effective  error  estimators 
exhibit  increasing  accuracy  whenever  a  proper  adaptive  approach  is  used. 

The  mentioned  flexibility  of  our  finite-element  formulation  for  the 
model  problem  incorporates  both  the  H-version  and  the  P-version  of  the  method. 
In  other  words,  we  may  fix  the  degree  of  the  elements  and  achieve  the  desired 
accuracy  by  refining  the  mesh,  or,  alternately,  we  may  fix  the  mesh  and 
increase  the  degree  of  the  elements  to  meet  the  accuracy  requirement.  For 
some  details  and  comparisons  of  these  two  versions  see  [14],  [15],  [41]. 

In  this  paper  we  shall  concentrate  on  the  H-version  for  which  the  approaches 
presented  here  are  closely  related  to  those  analyzed  for  the  linear  case  in 
[4  ],  [5  ],  [6  ],  [9  ],  [13]. 

We  begin  with  a  summary  of  the  definition  of  the  error-estimators  for 
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the  case  d  =  1  of  C°-element  of  degree  p.  The  corresponding  proofs 
win  be  given  elsewhere. 

For  given  parameter  values  x,  v  we  denote  by  u(s)  =  u(x,v)(s) 
the  exact  solution  of  our  problem.  Moreover,  for  a  specified  mesh  a 
of  (4.2)  and  fixed  degree  p  let  u(s)  =  u(A,p,X,v)(s) ,  O^s^  1,  be 
the  computed  finite  element  solution.  Then  the  error  estimator  is  obtained 
as  follows; 


(i)  On  every  sub-interval  I.,  i  =  l,...,n,  compute  the 
finite  element  approximation 

p+1 

(7.1)  w.(s)  =  I  z..  4».(t.(s)),  Vs  e  I., 

'  1=0  ’J  J  ^ 


of  the  original  problem  subject  to  the  boundary  conditions 


(7.2) 


(7.3) 


w^.(s^)  =  u(a,p,X,v)(s^),  £  =  i-l,i. 

(Here  the  notation  of  section  3  is  used).  Then,  for  each 
i,  1  <  i  <  n,  define 


=  [ 


(w;.(s) 


ii'(s))^ds]^/^ 


^i 


(ii)  On  every  sub-interval  compute  the  residual 

(7.4)  r.(s)  =  a(u'(s))'  +  b(u(s),x,v) 


and  the  quantity 
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(7.5)  p.  =  c  hP[|  r.(s)^ds]^/^ 


where  c  is  a  fixed  constant  (see  [13]).  For  small 
h.  we  have  p.  «  a.  and  hence  this  step  (ii)  can 
often  be  omitted. 


(iii)  With  the  quantities  (7.3)  and  (7.5)  the  error  indicator 
of  the  ith  sub-interval  is  defined  by 

o  0  1  /2 

(7.6)  n^(A)  =  [a^  +  pf]  ,  i  =  l,...,n. 

(iv)  Now  the  quantity 

n  2  1/2 

(7.7)  cU}=[l  nA^ri 

i=l  ^ 

is  the  desired  error  estimator  of  the  computed  solution. 

The  error  estimator  e  of  (7.7)  approximates  the  error 

(7.8)  llu(x,v)  -  u(a,p,x,v)11  ,  =  [  (u'(s)  -  u'(s))^ds]^'^^ 

H  (I)  Jq 

of  the  computed  finite  element  solution  in  the  H^-norm.  More  precisely, 
under  fairly  general  assumptions  it  is  possible  to  prove  that  when  the  exact 
solution  u(x,v)  exists  in  a  suitable  neighborhood  of  the  computed  solution 
u(a,p,x,v),  then 

llu(x,v)  -  u(a,p,x,v)11  ,  =  e(A)(l  +  o(l)),  as  e  ->•  0. 

H^I) 


(7.9) 


The  term  o(l)  depends  on  the  exact  solution.  As  we  saw  earlier  the 
existence  assumption  is  certainly  needed  if  we  want  an  estimate  of  the 
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form  (7.9).  More  generally,  suppose  that  we  trace  a  curve  v  =  =  const. 

Then,  for  a  >  1  it  can  be  shown  that  there  exists  an  exact  solution 
u(a,v)  such  that 


(7.10) 


u(a,P,X,v  )  '  u(x,vj|  I  1 
0  a 


X|“  =  e(A)(l  +  o(l))  as  e  ^  0. 


For  a  solution  curve  defined  by  v  =  the  situation  corresponding  to 
the  two  estimates  (7.9)  and  (7.10)was  schematically  indicated  in  Figure  6, 
The  error  estimator  (7.7)  is  composed  of  the  error  indicators  (7.6) 
for  the  individual  sub- intervals  of  the  given  mesh.  For  the  linear  case 
we  have  shown  in  [9  ]  that  an  (asymptotically)  optimal  mesh  is  achieved  if 
we  equilibrate  the  error  indicators,  that  is,  if  we  construct  a  mesh  for 
which  these  indicators  are  equal.  More  specifically,  we  call  a  sequence 
of  meshes  equilibrated  if 


(7.11) 


max  n.(A.  ) 

<  K,  Vk, 
min  n^(A,^)  -  ’ 


with  a  constant  K  that  is  independent  of  the  particular  mesh  Aj^.  These 
concepts  provide  the  basis  for  an  adaptive  mesh-refinement  procedure  which 
generates  a  sequence  of  meshes  that  satisfy.  (7.11)  and  for  which  the  error 
indicators  are  approximately  equal,  (see  [8],  [10],  [11],  [12],  [46]).  As  noted 
before,  the  equilibration  of  the  meshes  is  expected  also  to  increase  the 
effectivity  of  the  error  estimator  (see  [12],  [13]). 

Suppose  again  that  we  trace  a  curve  with  "  const.  As  we 

observed  in  Sections  3  and  6  the  position  of  a  turning  or  limit  point  on 
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this  curve  frequently  has  considerable  practical  interest.  More  specifi¬ 
cally,  we  are  interested  in  estimating  the  error  |Xy  -  x.|.|  indicated 
in  Figure  6.  Let  u^  =  u(a ,p,X.|.,v^)  be  the  computed  solution  at  the 
computed  critical  point.  Then  an  a  posteriori  estimate  of  Uj  - 
be  obtained  by  the  following  procedure: 

(i)  Compute  the  indicators  a.,  i  =  l,...,n,  of  (7.3) 
and  the  quantity  ’ 

(7.12)  a  =  [  I 

i=l  ^ 

(ii)  Let  (p  be  the  function  corresponding  to  the  tangent 
vector  (5.12)  at  the  computed  critical  point  and  form 

(7.13)  <  =  0.05 

(iii)  On  every  sub- interval  I.,  i  =  l,...,n,  compute  the 
finite  element  approximation  w.  of  the  form  (7.1) 
subject  to  the  boundary  conditiins 

(7.14)  ^  <{>(s^),  i  =  i-l,i. 


(7.16) 
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where  the  last  quantity  represents  a  linearization 
of  the  non-linear  function  of  x 


(7.17) 


[b(u.j.,X,VQ) 


b(iY,Xy,v^) Jk  $  ds. 


(v)  Then  the  quotient 


(7.18)  6  =  |i| 
is  the  desired  estimator. 

Analogously  as  in  the  case  of  the  estimator  of  the  solution  it  is 
possible  to  show  that 

(7.19)  |Xy  -  Xjl  =  e(l  +  o(l)),  as  E  0 
where  the  term  c(l)  depends  on  the  exact  solution. 


i 
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8.  Numerical  Experiments 

The  theory  outlined  in  the  previous  section  raises  several  questions, 
as,  for  instance,  the  following; 

(i)  In  practice,  how  effective  are  the  error  estimators  for 
the  computed  solution  and  the  critical  values,  especially 
in  view  of  the  fact  that  the  theory  of  these  estimators  is 
only  asymptotic  in  nature? 

(ii)  Is  the  quality  of  the  estimators  influenced  by  the  type  of 
nonlinearity,  as  exemplified  for  instance  by  the  three 
problems  (2.14)? 

(iii)  How  effective  is  the  adaptive  procedure  based  on  the  equili¬ 
bration  of  the  error  indicators,  and,  in  particular,  its  use 
in  combination  with  the  continuation  process  of  Section  5? 

In  order  to  provide  some  answers  to  these  and  related  questions  we  pre¬ 
sent  here  some  computational  results  obtained  with  our  model  problem. 

We  consider  first  the  strongly  nonlinear  rod  (2.14)  (i)  and,  in  particular, 
the  general  shape  of  its  equilibrium  surface.  Figure  9  correspoids  in 
essence  to  Figure  4  for  the  simple  framework  of  Section  3.  More  specifically, 
uniform  meshes  with  eight  subintervals  and  cubic  C°-elements  were  used,  and 
Figure  9  shows  several  solution  curves  for  v  =  const,  plotted  in  a  coordinate 
system  with  x  as  the  abscissa  and  the  norm 

(8.1)  11^)1,  =  niax  lu(x)l 

as  ordinate.  Note  that  in  contrast  to  the  case  of  the  framework  of  Figure  2 
the  solutions  now  have  two  turning  points.  The  first  bifurcation  point  occurs 
approximately  at  x  *  2.8954,  v  =  0,  and  the  second  one  at  x  *  6.5691, 
u  =  0.  For  the  linearized  equation  the  corresponding  first  two  eigenvalues 
are  '•  -  2.8954  and  \  -  6.1448. 
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Figure  10  gives  a  corresponding  picture  of  the  solution  curves  with 
the  moment  balance  M(0)  -  M(l]  (see  (2.16))  as  ordinate  in  place  of  (8.1). 
Figures  11  and  12  show  the  shape  of  the  rod  for  v  =  0.1  and  variable  x 
and  for  x  =  3.0  and  variable  v,  respectively. 

In  order  to  assess  the  effect! vity  of  the  error  estimation,  we  used 
as  "nearly  exact"  solution  the  computed  solution  with  a  uniform  mesh  of  10 
subintervals  and  quintic  C°-elements.  We  consider  the  curve  v  =  0.1  and 
on  it  the  five  target  points  X  =  1.0,  2.5,  3.5,  2.5,  1.6.  Table  1  gives 
the  computed  moment  balances  M(0)  -  M(l)  at  these  five  points  for  a 
uniform  mesh  with  eight  subintervals  and  C°-elements  of  different  order 


1 

3 

5 

1.0 

0.1203 

0.1206 

0.1206 

2.5 

0.5362 

0.5558 

1 

0.5558 

3.5 

2.595 

2.636 

2.636 

2.5 

2.787 

2.500 

2.500 

1.6 

2.845 

2.830 

2.830 

1-  -  -  -  - 

Table  1 

Table  2  shows  the  computed  error  estimators  and  the  "nearly  exact"  error 
norm  in  percent  of  the  norm  of  the  solution.  In  these  and  all  other  compu¬ 
tations  in  this  section  the  (small)  contribution  of  the  p^-terms  of  (7.5) 
was  neglected.  The  table  shows  clearly  the  high  quality  of  the  estimators. 

In  the  cases  p  =  4,  x  =  1.0  and  x  =  2.5  left  blank  in  the  table  the 
error  is  essentially  zero  within  machine  accuracy. 
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1 

1 

2 

3 

4 

Error  estim. 

1.0 

12.15% 

0.1137% 

0.004089% 

Error  norm 

12.13% 

0.1138% 

0.004089% 

Error  estim. 

2.5 

11.55% 

0.4306% 

0.01552% 

Error  norm 

12.04% 

1 

0.4310% 

0.01552% 

Error  estim. 

1 

3.5 

21 .58% 

1.597% 

1.123% 

0.5011% 

Error  norm 

1 

1 

23.91% 

1.900% 

1.391% 

0.5389% 

1 

Error  estim. 

2.5 

21.53% 

1.708% 

o 

o 

0.5013% 

Error  norm 

21 .88% 

1.865% 

1.285% 

0.5455% 

Error  estim. 

1 .6 

19.45% 

0.7632% 

0.8123% 

0.3171% 

Error  norm 

19.14% 

1.135% 

0.9576% 

0.3312% 

Table 

As  Figure  10  shows  there  are  two 
V  =  0.1.  With  a  mesh  of  eight  equal, 
occur  at 

(8.2)  \j  i  3.6546, 

'l 

In  order  to  assess  the  effect! vity  of 
value  we  used  for  the  comparison  five 
respectively. 


2 

turning  points  along  the  curve  for 
cubic  C°-elements  these  turning  points 

Xt  =  1.46749. 

*2 

our  error-estimators  of  these  critical 
and  ten  equal,  quadratic  C°-eleiTients , 
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Table  3  gives  the  result  of  these  computations.  Again  the  error  estimators 
perform  very  well.  The  table  also  shows  that  convergence  of  the  critical 
x-values  is  faster  than  that  of  the  corresponding  solution  errors.  In 
fact,  asymptotically  it  is  expected  that  the  critical  values  converge  with 
the  square  of  the  solution  errors  in  the  H^-norm. 


Limit 

Point 

No.  of 
Interv. 

Est.  of 
Solution 
Error 

Val  ue 

Limit  Point 

Error 

Estim. 

5 

13.11% 

3.699385 

0.04472 

0.05294 

10 

3.18% 

3.657982 

0.00317 

0.00399 

5 

7.205% 

1.470290 

0.002793 

0.002843 

'^2 

10 

1.61 3% 

1 .467897 

0.00040 

0.000357 

Table  3 


In  preparation  for  the  use  of  the  error  estimators  for  an  adaptive  con¬ 
struction  of  optimal  meshes,  we  show  in  Table  4  the  error  indicators  of  the 
solution  at  x  =  3.5,  v  =  0.1  for  uniform  C°-elements  of  varying  order. 


Evidently  for  lower  p  the  mesh  is  better  equilibrated  than  for 
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1 

2 

3 

4 

5 

1 

0.4085 

0.2849(-l) 

0.2525(-l) 

0.1133{-1) 

i  0.2058(-3) 

2 

0.1789 

0.2149(-1) 

0.2816(-2) 

0.3089(-3) 

.  0.2637(-4) 

3 

0.9545 

0.5290(-2) 

0.5770(-3) 

0.4769(-4) 

0.4212(-5) 

4 

0.7281 

0.1173{-2) 

0.2176(-3) 

0.7614{-5) 

0.9042{-6) 

5 

0.7281 

0.1173(-2) 

0.2176(-3) 

0.7614(-5) 

1 

,  0.9060(-6) 

6 

0.9545 

0.5290(-2) 

0.5770(-3) 

0.4769{-4) 

1  j 

0.4212(-5) 

7 

0.1789 

0.2149(-1) 

0.2816(-2) 

0.3089{-3) 

0.2637(-4) 

8 

1 

0.4085 

0.2849(-l) 

0.2525(-l) 

0.1133(-1) 

0.2058(-3) 

max 

min 

5.34 

24.3 

116. 

1 

i 

1 

1488. 

228. 

Table  4 

higher  values.  (In  the  case  p  =  5  in  Table  4  the  very  small  error  indicators 
for  the  intervals  with  index  i  =  4  and  i  =  5  are  distortet  by  round  off 
which  causes  their  values  to  be  too  large.)  The  results  of  Table  4  indicate 
that  an  adaptive  procedure  for  the  construction  of  the  mesh  should  be 
effective. 

As  Table  2  shows  the  error  may  vary  considerably  along  any  one  solution 
curve.  Suppose  that  we  wish  to  compute  the  curve  for  v  =  0.1  with  an 
accuracy  of  2Z,  For  this  we  begin  with  a  mesh  that  meets  this  requirement; 
in  our  case,  a  uniform  mesh  of  two  quadratic  C°-elements.  Then  during  the 


Table  5 
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In  order  to  compare  the  results  of  Tables  5  and  6  we  define  the 
work  per  continuation  step  as 

(8.3)  f^m  +  f2m^'^  units, 

where  f^  and  f2  are  the  number  of  function  and  derivative  evaluations, 
respectively,  and  m  is  the  number  of  degrees  of  freedom  (see  Section  4). 

For  the  computation  with  the  adaptive  procedure,  as  shown  in  Table  5, 
the  total  work  was  45,672  units  while  for  the  run  with  the  uniform  mesh, 
given  in  Table  6,  it  was  143,383  units.  Thus,  the  adaptive  process 
requires  only  31.9%  of  the  work  for  the  fixed  uniform  mesh.  At  the  same 
time,  a  maximum  of  eight  elements  were  needed  during  the  adaptive  process 
while  the  fixed  uniform  mesh  with  the  same  accuracy  required  twelve  elements. 
These  results  are  typical  for  problems  of  this  type.  The  error  usually 
varies  considerably  along  any  solution  curve  and  it  is  practically  impossible 
to  design  a  fixed  uniform,  or,  for  that  matter,  non-uniform  mesh  for  which 
it  can  be  guaranteed  that  the  error  remains  below  a  prescribed  accuracy. 

Only  an  adaptive  process  which  modifies  the  mesh  in  response  to  the  computed 
error  estimates  can  accomplish  this  task.  In  our  case  we  used  only  a  mesh- 
refinement  strategy,  but,  obviously,  it  is  also  possible  to  incorporate  de¬ 
refinement  when  the  estimates  fall  below  a  lower  threshold. 

It  may  be  of  interest  to  analyze  why  the  fixed  uniform  meshes  do  not 
perform  well  in  the  higher  parts  of  our  solution  curve  for  v  =  0.1.  For 
this  we  use  a  uniform  mesh  of  eight,  cubic  C°-elements  and  give  in  Table  7a 
the  error  indicators  for  the  first  three  intervals  at  various  target  points 
along  the  curve  v  =  0.1.  In  addition.  Table  7b  shows  the  actual  values  of 
the  function  u  and  its  derivative  u'  at  the  corresponding  nodal  points 
and  the  midpoints  of  these  intervals.  Clearly,  in  the  upper  parts  of  the 
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curve  boundary  layers  develop  near  the  ends  of  the  rod  and  this  is  re¬ 
flected  by  the  error  indicators. 


2.5 

!  3.5 

2.5 

X  =  1.4675 
Limit  Point 

1 

0.1703(-4) 

1 

0.2524(-l) 

0.3200(-l) 

0.2778(-l) 

2 

0.2358(-4) 

0.2816(-2) 

0.3924(-2) 

0.5590(-2) 

1 

3 

0.2218(-4) 

0.5770(-3) 

0.7445(-3) 

0.2720(-2) 

Table  7a 


X  =  2.5  X  =  3.5  X  =  2.5  x  =  1.4675 


X 

u 

u' 

u 

u' 

u 

u' 

u 

u* 

0 

0 

0.5707 

0 

2.730 

0 

11.17 

0 

20.74 

0.0625 

0.03482 

0.5409 

0.1648 

2.515 

0.5833 

7.639 

1.048 

13.20 

0.1250 

0.06726 

0.4946 

0.3106 

2.133 

0.9728 

5.013 

1.702 

8.126 

0.1875 

0.09632 

0.4335 

0.4304 

1 .703 

1.241 

3.644 

2.091 

4.597 

0.2500 

0.1212 

0.3606 

0.5240 

1.299 

1  .438 

2.714 

1 

2.311 

2.636 

0.3125 

0.1412 

0.2782 

0.5937 

0.9358 

1.585 

1.989 

2.442 

1.628 

0.3750 

0.1558 

0.1892 

0.6417 

0.6044 

1 .687 

1 .297 

2.521 

0.9527 

Table  7b 


63 


So  far  we  discussed  only  computations  with  C°-elements.  In  order 
to  assess  the  influence  on  the  error  of  different  degrees  of  smoothness 
we  selected  three  cases  with  the  same  number  of  degrees  of  freedom  (4.5), 
namely  m  =  10.  The  corresponding  results  for  the  solution  curve  defined 
by  V  =  0.1  in  the  range  0  _<  x  _<  1.467  up  to  the  second  limit  point 
are  given  in  Table  8: 


Smoothness 

Degree 

No.  of 
Elements 

Max.  Error 
in  Range 

Work 

Units 

C° 

1 

11 

23.81% 

22,243 

3 

5 

3.939% 

20,236 

C2 

5 

1 

3 

0.9976% 

20,919 

Table  8 


The  effectivity  of  the  P-version  inherent  in  these  results  corresponds  to  that 
in  the  linear  case  discussed  in  [2],  [15]. 

All  results  given  so  far  in  this  section  concerned  the  strongly  non¬ 
linear  rod  (2.14)  (i).  The  results  are  essentially  analogous  to  that  for 
the  mildly  non-lirear  case.  However,  the  turning  points  now  occur  for 
larger  values  of  x.  This  can  be  seen  from  Figure  13  where  we  give  the  first 
bifurcating  solution  for  v  =  0.0  in  the  three  cases  (2.14).  Ten  equal,  quintic 
C°-elements  were  used  in  each  case.  The  effectivity  of  the  error  estimators 


also  remains  the  same  in  the  two  cases  a  comparison  of  Tables  2  and  9  shows. 


X 


1 


2 


3 


4 


Error  estim. 

Error  norm 

3.0 

11.45% 

11 .61% 

.5184% 

.5187% 

.01964% 

.01965% 

.6520(-3) 

.6535(-3) 

1 

1 

Error  estim. 

5.0 

11.87% 

.7603% 

.1023% 

j 

.3493(-2) 

Error  norm  | 

n.92% 

. 7672% 

I  .1024% 

i 

.2906(-2) 

Error  estim. 

1 

7.0 

12.73% 

.8964% 

1  .2426% 

1 

.3202(-l) 

Error  norm 

i 

12.83% 

.9296%  ' 

.2447%  1 

.3202(-l) 

Error  estim. 

5.0 

14.49% 

1.790% 

.1260% 

.8446(-l) 

Error  norm 

1 

1  12.26% 

i 

1.797% 

.1518% 

.8448(-l) 

Error  estim. 

3.0 

13.67% 

1 

1.166% 

.2288% 

.1595(-1) 

Error  norm 

1 

12.81% 

1.189% 

.2296% 

.1850(-1) 

Table  9 


In  suimary  the  computational  results  of  this  section  provide  support 
for  the  following  claims: 


a)  The  a  posteriori  estimators  of  the  errors  of  the  solution  and 
the  critical  values  are  very  reliable  for  different  types  of 
nonlinearities  and  all  degrees  of  elements. 

b)  Higher  order  elements  perform  in  general  better  than  the  lower 
order  ones . 

c)  The  adaptive  procedures  for  modifying  the  mesh  during  the  con¬ 
tinuation  process  is  very  effective  and  significantly  increases 
the  effectivity  of  the  computations. 


i 


9.  Spurious  Solutions 

As  mentioned  in  the  Introduction,  for  nonlinear  problems  of  the  type 
considered  here  the  equilibrium  surface  of  the  discretized  equations  often 
has  a  different  number  of  connected  components  than  that  of  the  original 
equations.  The  components  which  do  not  approximate  reasonably  well  the 
exact  solutions  in  any  of  their  parts  have  been  called  spurious,  or 
numerically  irrelevant  solutions.  They  have  been  observed  by  many  authors 
(see  eg.  [20],  [28]  and  the  references  given  there). 


The  phenomenon  is  easily  observed  in  the  case  of  the  classical  Euler 
rod  (2.14)  (iii).  For  v  =  0  it  can  be  shown  that  for  x  <  0  there  exists 
no  other  solution  than  u  h  0.  This  is  physically  plausible  since  under 
pure  tension  the  rod  should  remain  straight.  In  fact.  Figure  13  showed 
that  the  first  non-zero  solution  bifurcates  from  the  trivial  solution  at 
X  =  4.940,  and  hence  that  we  should  expect  non-zero  solutions  for  v  =  0 
only  when  x  exceeds  that  value.  Nevertheless,  the  finite  element  method 
does  produce  non-zero  solutions  for  certain  negative  values  of  x.  Figure 
14  shows  this  for  the  cases  of  uniform  meshes  with  two,  four,  and  six 
linear  C°-elements,  respectively.  The  limit  points  with  respect  to  x  for 
these  three  paths  as  well  as  for  that  one  obtained  with  eight  linear  elements, 
are  given  in  Table  10. 


Number  of  \j 
Elements 

2  69.629 
4  47.404 
6  64.246 
8  110.907 


Table  10 
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In  this  simple  case  it  is  easy  to  distinguish  between  the  correct 
solutions  for  x  >  0  and  the  spurious  ones  for  x  <  0.  In  general, 
however,  the  situation  is  more  complicated  and  it  is  not  readily  clear 
how  to  distinguish  between  the  correct  and  the  spurious  solutions.  Under 
fairly  general  conditions  the  finite  element  solutions  converge  to  the 
solutions  of  the  infinite  dimensional  problem.  Hence,  in  general,  the 
spurious  solutions  have  to  disappear  toward  infinity  when  the  mesh  is  re¬ 
fined.  But  this  need  not  happen  monotonical ly  as  Table  10  shows. 

In  line  with  our  overall  philosophy,  we  are  interested  in  computing 
an  approximate  solution  for  which  the  a  posteriori  estimate  indicates  a 
deviation  of  at  most  10-15^  from  an  exact  solution.  In  other  words,  we  will 
accept  only  those  segments  of  any  computed  solution  path  for  which  the 
estimate  does  not  exceed,  say,  12%.  The  spurious  solutions  should  carry  a 
much  larger  error  and  hence  could  be  rejected  on  that  basis.  This  raises 
the  question  whether  our  error  estimators  will  succeed  in  identifying 
properly  such  large  errors  even  though  they  were  justified  theoretically 
only  asymptotically  for  small  errors. 

The  following  tables  show  that  this  is  indeed  the  case.  More  specifi¬ 
cally,  Table  11  gives  some  error  estimates  for  the  spurious  solution  path 
obtained  with  a  uniform  mesh  of  eight  linear  C°-elements. 


< 

\ 

Hull  1 
h' 

Error 

Estimate 

-128.502 

15.815 

26.9% 

-125.507 

15.688 

27.0% 

-122.513 

1 

15.553 

27.2% 

-119.522 

15.409 

27.5% 

-116.534 

15.250 

27.8% 

-113.554 

15.064 

28.4% 

-110.907* 

14.758 

30.4% 

-112.192 

14.632 

32.5% 

-115.796 

14.581 

35.0% 

-119.599 

14.585 

36.8% 

-123.491 

14.611 

38.4% 

-127.417 

14.651 

39.9% 

-131.359 

14.699 

41.3% 

*  Limit  Point 

Table  n 

Since  u  ;  0  is  the  only  exact  solution  for  negative  A  the  value  of 

I  lull  represents  here  the  actual  error  in  the  H^-norm.  In  other  words, 
h' 

the  actual  error  is  100%,  and,  not  surprisingly,  the  error  estimators  only 
have  an  efficiency  index  of  about  0.27-0.41.  Yet  in  line  with  the  earlier 
discussion  any  solution  path  such  as  this  one  with  errors  exceeding  the 
12%  threshold  should  be  rejected. 
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On  the  other  hand,  Table  12  gives  the  error  estimates  for  the  "correct" 
solution  path  bifurcating  from  the  positive  X-axis  at  about  x  =  5. 

Again  a  uniform  mesh  with  eight  linear  elements  was  used. 


Error 

X  Estimates 


5.1390 
5.3264 
5.5972 
5.9553 
6.3551 
6.7940 
7.3081 
7.8960 
8.5529 
9.2715 
10.042 
10.856 
11 .704 
12.579 


11.43% 
11.43% 
11.41% 
11.41% 
11.41% 
11.42% 
11.45% 
11.49% 
11.54% 
11 .62% 
11.71% 
11 .82% 
11.95% 
12.09% 


Table  12 


If  we  allow  a  maximum  error  of  12%  then  the  entire  path  segment  up  to 
\  =  12.579  should  be  accepted.  On  the  other  hand,  if  the  maximum  error 
is  10%  then  this  path  must  be  rejected  and  higher  order  elements  are  needed. 
Hence,  the  question  of  accepting  or  rejecting  a  path  segment  rests  solely  on 
our  error  estimate. 


71 


Part  II:  Outlook  on  More  General  Linear  Problems 

The  results  discussed  in  Part  I  above  are  as  yet  still  rather  re¬ 
strictive  in  their  applicability.  Clearly,  we  should  hope  to  be  able  to 
extend  them  to  boundary  value  problems  in  several  space  dimensions  of  the 
type  occurring  in  nonlinear  structural  analysis.  But  there  are  other 
possible  extensions  as  well.  In  this  part  we  illustrate  some  of  the  desired 
extension  by  means  of  results  already  obtained  for  the  linear  case. 

More  specifically,  we  begin  with  some  computational  examples  of  the 
use  of  our  adaptive  finite  element  solver  FEARS.  This  system  is  based  on 
our  results  cited  in  the  Introduction  about  error  estimators  and  adaptive 
procedures  for  two-dimensional  elliptic  problems.  Then  we  turn  to  a  simple 
parabolic  model  problem  to  show  the  possibility  of  applying  these  ideas 
also  in  the  case  of  certain  time-dependent  problems. 

As  indicated  before  the  material  in  this  Part  II  is  strictly  intended 
only  to  indicate  the  type  of  results  we  would  hope  to  extend  also  to  the 
nonlinear  case.  For  more  details  on  these  linear  problems  we  refer  to  the 


cited  articles. 
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10.  An  Adaptive  Finite  Element  Solver 

For  linear  elliptic  problems  in  two-space  dimensions  the  theory  of 
a  posteriori  error  estimates  has  been  advanced  to  a  point  (see  eg.  [4], 

[5],  [11])  that  it  became  possible  to  implement  a  general  finite  element 
solver  which  incorporates  the  ideas  outlined  here,  (see  eg.  [8],  [36],  [46]). 
More  specifically,  the  FEARS  system  [Finite  Element  Adaptive  Research  Solver] 
was  developed  at  the  University  of  Maryland  and  is  available  in  FORTRAN  for 
use  on  Uni  vac  1100  series  computers  and  on  the  CDC-7600.  It  has  the  following 
principal  features: 

(i)  FEARS  constitutes  an  appl ications- independent  finite  element 
solver  for  a  certain  class  of  two-dimensional,  linear,  elliptic 
boundary  value  problems  defined  by  a  weak  mathematical  formu¬ 
lation.  At  present  curvilinear  elements  of  degree  1  are  used. 

(10.1)  (ii)  Adaptive  approaches  are  employed  extensively.  The  a  posteriori 
estimates  developed  by  us  (loc.cit)  are  used  to  control  the 
adaptive  processes  and  to  provide  a  solution  with  a  near 
optimal  error. 

(iii)  In  the  systems  design,  advantage  was  taken  of  the  natural 
modularity  of  the  finite  element  method  (see  [45],  [46]). 

Some  results  obtained  with  the  system  were  described  in  [12].  We  pre¬ 
sent  here  some  further  results  which  illustrate  the  effect! vity  of  the 
adaptive  approach. 

Example  1 :  We  consider  the  elasticity  problem  for  an  isotropic  homogeneous 
material  with  Young's  modulus  E  =  1  and  Poisson's  ratio  v  =  0.3  on  the 
domain  shown  in  Figure  15  where  also  the  boundary  conditions  are  specified. 

The  energy  norm  is  used  for  the  error  estimation.  The  basic  mesh  consisted 
of  eight  equal,  bilinear  element  and  Figure  16  a-j  provides  a  sequence  of 
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nine  consecutive  meshed  generated  by  FEARS.  The  behavior  of  the  (energy- 
norm)  error  for  this  problem  in  dependence  of  the  degrees  of  freedom  N 
is  shown  in  Figure  17.  Note  that  N”^^^  is  the  fastest  possible  rate 
of  convergence  while  the  rate  for  uniform  meshes  is  only  because 

of  the  presence  of  the  singularity.  Hence,  we  see  from  Figure  17  that  the 
adaptive  procedure  provides  us  with  a  sequence  of  solutions  with  the  fastest 
possible  rate  of  convergence.  Finally  Figure  18  gives  the  effectivity  index 
(1.1)  for  the  estimator  used  in  the  process.  Obviously,  for  errors  of  10% 
or  better  the  accuracy  is  very  high. 

Example  2:  As  a  second  elasticity  problem  we  consider  the  behavior  of 
isotropic,  homogeneous  material  with  E  =  1,  v  =  0  on  a  ring-domain 
defined  by  concentric  rings  with  radii  0.2  and  1,  respectively.  On  the 
inner  ring  we  specify  the  boundary  conditions  u  =  1,  v  =  0  and  on  the 
outer  ring  u  =  0,  v  =  0.  Because  of  the  symmetries  of  the  solution, 
we  restrict  the  computation  to  a  quarter  of  the  ring  with  obvious  boundary 
conditions  on  the  radial  part  of  the  boundary.  As  in  Example  1  the  error 
is  measured  in  the  energy  norm.  The  sequence  of  meshes  generated  by  FEARS 
is  shown  in  Figure  19  and  the  corresponding  error  behavior  is  given  in  Table 
13. 


ENERGY  NORM  OF  THE  ERROR  IN  % 


Figure  17 


ADAPTIVE  S 
MESH  y\ 
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No.  of 
Elem. 

Degrees 

of 

Freedom 

j 

Actual 

Error 

llellE 

Error 

Estim. 

1 

Rel .  Error 

100||e||/||u|| 

Effect. 

6 

4 

4 

0.2495 

0.2210 

29.9% 

0.885 

16 

24 

0.1392 

0.1135 

1  16.7% 

0.815 

48 

80 

0.0784 

0.0738 

9.4% 

0.942 

145 

225 

0.0439 

0.0434 

5.3% 

0.988 

Table  13 

Once  again  we  see  that  the  effectivity  index  0  of  (1.1)  has  practically 
acceptable  values  when  the  accuracy  of  the  solution  is  in  the  range  of  5-10%. 


11.  A  Parabolic  Problem 


The  ideas  and  approaches  discussed  here  can  also  be  applied  to  time- 
dependent  parabolic  problems.  We  restrict  our  comments  to  the  simple 
model  problem 


(11.1) 


au  _  ,  a  u 


0<t<T,  0<x<l 


with  the  initial  conditions 


(11.2)  u(0,x)  =  g(x) 

and  boundary  conditions 


(11.3) 


U(t,0)  =  U(t,l)  a  0. 


The  computations  involved  in  the  method  of  lines  consist  of  the  following 
two  parts: 

(i)  Use  a  finite-element  discretization  of  the  space  variable 
to  obtain  an  initial  value  problem  for  a  system  of  ordinary 
^11  differential  equations. 

(ii)  Solve  this  system  numerically  by  means  of  some  adaptive  solver. 


The  discretizations  in  space  and  time  have  a  different  character.  How¬ 
ever,  the  objective  is  to  equilibrate  these  two  errors.  We  sketch  here  the 
principal  idea  of  the  error  estimate  of  the  finite  element  discretization 
and  refer  for  more  details  to  [ 3  ]. 

If  the  time-derivative  du/dt  on  the  left  of  (11.1)  were  given,  then 
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the  problem  becomes  elliptic  and  we  know  how  to  estimate  the  error  in  the 
energy  norm 

(11.5)  ||e||2  =  a(||)2  dx. 

It  is  possible  to  show  that  with  the  values  of  du/dt  defined  by  an  approximate 
solution  of  the  system  of  differential  equations  (11.4)  (ii)  we  still 
obtain  in  this  way  a  correct  error  estimator  with  an  effectivity  index  that 
converges  to  one  as  the  error  tends  to  zero  (see  [3]). 

We  illustrate  this  approach  on  our  model  problem  (11 .1 )-(ll .3)  with 
g(x)  =  sin  ttx,  a  =  0.3  and  T  =  1  for  which  the  differential  equations 
(11.4)  (ii)  can  be  solved  exactly.  Table  14  shows  the  values  of  the  effecti¬ 
vity  index  in  the  energy  norm  (11.5)  at  the  time  steps  t  =  1/3,  2/3,  1  in 
dependence  of  the  number  n  of  linear  elements  in  a  uniform  mesh  on  the 
X- interval  [0,1]. 


1/3 

2/3 

1 

4 

0.911258 

0.814769 

0.713659 

8 

0.977357 

1 

!  0.948020 

1 

0.910430 

16 

0.994310 

0.986583 

0.976027 

32 

0.998575 

0.996618 

0.993897 

64 

0.999643 

0.999152 

0.999614 

128 

0.999909 

0.999786 

0.999892 

Table  14 
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Evidently  the  effect! vity  index  converges  to  one  as  n  tends  to  infinity. 

The  above  results  assume  that  the  differential  equations  (11.4)  (ii) 
are  solved  explicitly.  In  practice,  of  course,  this  is  not  the  case  and 
we  have  to  use  one  of  the  modern  ODE-solvers  as,  for  instance,  LSODI 
which  handles  systems  in  the  linearly  implicit  form 

(11.6)  A(t,y)y  =  b(t,y). 

The  solver  selects  adaptively  the  steps! ze  and  order  so  as  to  meet  a 
specified  tolerance  t.  Hence,  our  objective  will  be  to  choose  the 
tolerance  x  as  a  function  of  the  space  discretization  and  the  step  and 
order  selection  criterion  of  the  ODE  solver  such  that 

(i)  the  error  of  the  semi  discrete  method  is  preserved 
at  the  least  cost,  and 

(ii)  the  effect! vity  of  the  estimators  is  maintained. 

These  two  objectives  can  be  achieved  with  considerable  success.  We  refer 
for  details  to  [ 3  ] . 
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